<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Sparse covariance estimation for Gaussian variables</title>
<link rel="canonical" href="/Users/mcgrant/Projects/CVX/examples/log_exp/html/sparse_covariance_est.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Sparse covariance estimation for Gaussian variables</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Jo&Atilde;&laquo;lle Skaf - 04/24/08</span>
<span class="comment">% (a figure is generated)</span>
<span class="comment">%</span>
<span class="comment">% Suppose y \in\reals^n is a Gaussian random variable with zero mean and</span>
<span class="comment">% covariance matrix R = \Expect(yy^T), with sparse inverse S = R^{-1}</span>
<span class="comment">% (S_ij = 0 means that y_i and y_j are conditionally independent).</span>
<span class="comment">% We want to estimate the covariance matrix R based on N independent</span>
<span class="comment">% samples y1,...,yN drawn from the distribution, and using prior knowledge</span>
<span class="comment">% that S is sparse</span>
<span class="comment">% A good heuristic for estimating R is to solve the problem</span>
<span class="comment">%           maximize    logdet(S) - tr(SY)</span>
<span class="comment">%           subject to  sum(sum(abs(S))) &lt;= alpha</span>
<span class="comment">%                       S &gt;= 0</span>
<span class="comment">% where Y is the sample covariance of y1,...,yN, and alpha is a sparsity</span>
<span class="comment">% parameter to be chosen or tuned.</span>

<span class="comment">% Input data</span>
rand(<span class="string">'state'</span>,0);
randn(<span class="string">'state'</span>,0);
n = 10;
N = 100;
Strue = sprandsym(n,0.5,0.01,1);
R = inv(full(Strue));
y_sample = sqrtm(R)*randn(n,N);
Y = cov(y_sample');
alpha = 50;

<span class="comment">% Computing sparse estimate of R^{-1}</span>
cvx_begin <span class="string">sdp</span>
    variable <span class="string">S(n,n)</span> <span class="string">symmetric</span>
    maximize <span class="string">log_det(S)</span> <span class="string">-</span> <span class="string">trace(S*Y)</span>
    sum(sum(abs(S))) &lt;= alpha
    S &gt;= 0
cvx_end
R_hat = inv(S);

S(find(S&lt;1e-4)) = 0;
figure;
subplot(121);
spy(Strue);
title(<span class="string">'Inverse of true covariance matrix'</span>)
subplot(122);
spy(S)
title(<span class="string">'Inverse of estimated covariance matrix'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
 
Successive approximation method to be employed.
   For improved efficiency, SDPT3 is solving the dual problem.
   SDPT3 will be called several times to refine the solution.
   Original size: 503 variables, 223 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  1/  1 | 3.471e+00  6.978e-01  0.000e+00 | Solved
  1/  1 | 3.433e-01  7.438e-03  0.000e+00 | Solved
  1/  1 | 4.190e-03  1.097e-06  0.000e+00 | Solved
  0/  0 | 0.000e+00  0.000e+00  0.000e+00 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -31.2401
 
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="sparse_covariance_est__01.png" alt=""> 
</div>
</div>
</body>
</html>